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SECURITY  CLASSIFICATION  OF  THIS  PACE 


Time  domain  waveform  modeling  studies  of  short  period  and  high  frequency 
have  been  carried  out  in  an  attempt  to  improve  event  discrimination  capabilities 
with  regional  data.  The  purpose  was  to  develop  discriminants  with  a  well 
understood  physical  basis  which  could  be  reliably  transported  to  different  areas 
of  the  world.  The  study  of  P^  involved  adapting  methods  developed  to  fit  long 
period  data  for  observations  with  periods  as  short  as  two  seconds.  To  accomplish 
this,  we  have  modified  the  wavenumber  integration  method  for  efficient  computation 
of  Pnl.  Ix>ng  period  Pn]>)  is  generally  computed  assuming  a  simple  layer  over  half 
space  modgl^^ Numerical  tests  were  performed  to  determine  how  the  model  needs 
to  be  generalized  to  fit  much  shorter  period  data-.  It  was  found  thatr^velocity 
gradients  like  those  typically  observed  near  the  free  surface  were  of  great 
importance.  Gradients  at  lower  depths  such  as  at  the  crust-mantle  transition 
are  of  much  less  significance.  A  simple  layered  model  which  accurately  predicts 
regional  P^;  propagation  from  NTS  to  the  surrounding  digital  stations  was  developed. 
It  deterministically  explains  the  waveshape  of  the  Pff  and  approximately  the  first 
15  seconds  of  high  amplitude  energy  reflecting  past  critical  angle  from  the 
mantle.  The  modeling  study  of  short  period  Pjf  has  led  to  the  development  of  an 
effective  and  physically  based  regional  discriminant;  the  underlying  purpose  of 
this  project.  It  has  been  found  that  the  waveshape  of  the  Pn  is  stable,  similar 
at  most  stations  and  straightforward  to  model.  Simple  correlation  of  observed 
Pn  waveforms  with  the  expected  average  waveform\or  explosions  is  sufficient  to 
discriminate  events.  In  preceding  reports,  it  wa£  established  that  high  frequency 
Pn  appeared  more  to  be  a  turning  ray  in  the  mantle  rather  than  a  classical  head 
wave.  It  was  also  found  that  pPn  arrived  anomalously  late.  In  the  latest  work, 
development  of  Pn  through  its  crossover  was  modeled,  resulting  in  a  structure 
appropriate  for  paths  from  NTS  to  Las  Vegas.  It  was  found  that  propagation  of 
P  to  regional  distances  is  best  modeled  in  terms  of  smooth  continuous  positive 
gradients  with  no  first  order  discontinuities.  A  model  for  the  anomalous  pPn 
behavior  was  developed  in  terms  of  spall. 
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This  is  the  final  report  on  a  two  year  effort  to  develop  effective  dis¬ 
criminants  between  explosion  and  earthquake  seismic  sources  which  make  use  of 
regional  seismic  recordings.  Regional  body  phases  have  proved  to  be  among  the 
most  difficult  to  analyze  or  model  in  a  deterministic  fashion  in  the  field  of 
seismology.  At  high  frequencies,  the  effects  of  scattering  in  the  crust  become 
so  intense  that  only  statistical  properties  of  waveforms  are  meaningful.  We 
have  undertaken  two  investigations  in  this  work  directed  at  overcoming  these 
difficulties  insofar  as  they  relate  to  seismic  discrimination.  The  first  was 
an  analysis  of  the  composition  of  the  portion  of  the  regional  signal  named  Pni 
by  Helmberger  and  Engen  (1980).  It  consists  of  approximately  the  first  60  to 
90  seconds  of  compressional  signal.  In  the  long  period  channel  this  includes 
the  long  period  head  wave  and  the  PL  wave  which  accounts  for  the  choice  of 
nomenclature.  In  the  short  period  channel,  this  window  would  encompass  high 
frequency  Pn  and  Pg.  The  second  study  is  an  attempt  to  develop  a  discriminant 
based  on  the  waveform  of  high  frequency  Pn.  These  two  studies  have  several 
features  in  common.  The  first  is  that  they  made  use  of  essentially  the  same 
data  base;  that  is,  records  from  the  regional  digital  seismic  network  recently 
set  up  around  the  Nevada  Test  Site  by  a  number  of  different  institutions.  Another 
is  that  they  both  take  advantage  of  the  fact  that  some  properties  of  regional 
signals  become  stable  and  deterministic  when  several  groups  of  signals  are 
appropriately  averaged  together.  This  technique  makes  the  best  use  of  those 
properties  of  regional  phases  that  are  meaningful  in  a  statistical  sense.  Finally, 
each  takes  advantage  of  an  exception  to  the  generality  that  regional  phases  are 
too  unstable  to  model  with  synthetic  seismograms.  This  is  apparently  not  true 
at  longer  period  or  at  very  short  times  at  the  beginning  of  the  signal. 


Three  interim  reports  have  been  distributed  chronicling  the  progress  of 
these  investigations.  We  begin  with  a  brief  summary  of  previous  work  on  short 
period  Pnl  which  is  relevant  to  the  final  outcome  of  the  investigation.  We  then 
discuss  the  final  results  on  this  portion  of  the  study.  We  present  a  similar 
synopsis  of  previous  work  and  finally  results  on  studies  of  short  period  Pn 
waveforms . 
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initial  work  oh  short  perioj? 

The  initial  motivation  for  Investigating  the  physical  nature  of  the  short 
period  energy  in  Pnl  came  from  the  recognition  that  it  would  be  very  useful  to 
develop  ways  to  passively  calibrate  regional  arrays,  When  regional  networks  are 
installed  in  new  areas  for  monitoring  purposes,  it  can  be  anticipated  that  they 
will  be  of  limited  use  because  of  uncertainties  about  regional  wave  propagation 
in  the  place  where  they  are  located.  However,  they  will  immediately  begin  to 
record  small  to  moderate  sized  earthquakes  at  regional  distances.  Events  in  the 
magnitude  range  of  3  to  5  occur  regularly  in  every  continental  area  and  most  of 
these  events  will  be  reliably  identified  as  either  explosions  or  earthquakes  by 
standard  methods.  Figure  1  shows  a  typical  broad  band  regional  signal  from  an 
earthquake  recorded  at  a  modern  digital  station.  These  are  the  compressional 
signals  from  the  Idaho  earthquake  of  October  29,  1983  as  recorded  at  RSSD.  Even 
though  the  magnitude  was  only  5.4,  the  event  produced  a  clear  signal  800  km  away 
on  all  three  channels.  The  long  period  trace  is  very  simple  and  can  be  easily 
modeled  using  standard  long  period  methods.  The  intermediate  trace  is  much 
richer  in  character.  The  Pn  onset  is  weak  followed  by  two  stronger  arrivals 
which  are  the  depth  phases  pPn  and  sPn.  There  is  a  distinct  large  amplitude 
arrival  about  twenty  seconds  into  the  record  with  a  simple  enough  character  to 
warrant  an  attempt  at  deterministic  modeling.  An  examination  of  the  particle 
motion  shows  that  some  scattering  by  non  stratified  structure  is  occurring. 
The  bottom  short  period  channel  exhibits  a  typical  high  frequency  regional  P 
wave.  There  is  a  high  frequency  Pn  followed  by  Pg.  The  relationship  between 
Pg  onset  and  the  phase  on  the  intermediate  channel  is  intriguing.  The  particle 
motion  in  the  short  period  channel  indicates  strong  scattering  with  equal  energy 
appearing  on  the  tangential,  radial  and  vertical  traces.  The  clear  implication 
of  Figure  1  is  that  there  is  a  breakdown  of  deterministic  character  of  the  signal 
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Figure  1 .  Long,  intermediate  and  short-period  vertical  Channels,  P  waves  from  the  Idaho  earthquake 
of  October  29,  1983  (M-5.4)  as  recorded  at  RSSD  (R-800  Km).  Seismograms  are  displayed  on  the 
same  time  scale.  First  arrival  is  Pn.  The  onset  of  large  amplitude  waves  is  the  arrival  time  of  Pf 
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descending  down  the  page.  However,  it  would  appear  that  if  the  top  two  traces 
could  be  deterministically  modeled,  they  would  serve  as  a  useful  guide  to 
interpretation  of  the  short  period  signal. 

The  basis  of  what  we  suggested  in  our  first  report  (Burdick  and  Helmberger, 
1987)  was  that  if  broad  band  signals  recorded  at  newly  installed  regional  stations 
were  routinely  modeled,  a  predictive  capability  could  be  developed  for  given 
paths.  Synthetic  long  and  intermediate  period  seismograms  could  be  computed  for 
any  path  from  which  a  reference  event  was  available.  The  utility  of  these 
synthetics  is  that  they  help  to  identify  when  phases  high  in  discrimination 
potential  such  as  depth  phases  or  shear  phases  arrive.  Unidentified  seismic 
events  are  generally  small  and  recorded  only  on  the  short  period  channel.  Longer 
period  synthetics  could  be  computed  for  them  using  the  calibration  models  developed 
over  time  and  discrimination  capability  could  be  Increased. 

We  discussed  our  forward  modeling  efforts  on  the  data  in  Figure  1  in  our 
first  report.  Some  of  the  results  are  shown  in  Figure  2.  The  short  period 
synthetics  shown  were  computed  using  the  formalism  outlined  in  Helmberger  and 
Engen  (1980)  along  with  their  layer  over  a  half  space  model  and  the  correct 
intermediate  period  instrument  response.  The  synthetics  show  a  series  of 
resonating  high  amplitude  arrivals  which  are  also  present  to  some  extent  in  the 
data.  We  named  these  phases  crustal  resonance  phases  and  examined  their  composition 
in  terms  of  generalized  rays  As  might  have  been  suspected,  these  resonances 
are  associated  with  successively  higher  orders  of  reverberations  in  the  crust. 
Interestingly,  the  second  resonance  phase  is  not  associated  with  the  second  order 
of  reverberation,  but  with  rays  that  have  at  least  one  converted  S  leg.  Thus, 
a  comparison  of  the  first  to  the  second  resonance  might  have  discrimination 
potential.  To  compute  complete  Pnl  synthetics,  Helmberger  and  Engen  (1980) 
typically  included  4094  generalized  rays  in  the  ray  set.  By  examining  small 
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Figure  2.  A  comparison  between  observed  and  synthetic  intermediate  period  P  wave  seismograms 
for  both  vertical  and  radial  components  from  the  Idaho  earthquake  of  October  29,  1983  (M-5.4). 
Crustal  resonance  phases  are  identified.  The  synthetics  were  computed  using  the  methodoIORv  of 
Helmberger  and  Engen  (1980). 
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subsets  of  rays,  we  established  that  the  first  resonance  phase  is  dominated  by 
the  six  rays  P,  S,  PmP ,  pPmP,  SmS  and  sSmS .  The  S  rays  contribute  through  the 
phenomenon  of  S  to  P  diffraction.  The  second  resonance  is  controlled  by  the 
eight  rays  PmS ,  SmP,  pSmS ,  sPmP,  PmPSmS ,  SmSPmP,  PmPSmP  and  PmSPmP.  We  performed 
additional  studies  which  indicated  that  the  Pn  signals  behave  more  like  turning 
rays  than  true  head  waves,  but  this  is  more  relevant  to  the  second  of  our 
investigations  as  discussed  below. 

To  further  establish  the  existence  and  stability  of  resonance  phases  at 
even  higher  frequencies,  we  attempted  stacking  high  frequency  signals  from  NTS 
explosions  as  observed  at  the  DWWSSN  station  ALQ.  The  results  were  discussed 
in  our  second  report  (Burdick  et  al.  ,  1988).  An  example  of  the  results  is  shown 
in  Figure  3.  The  crustal  resonances  emerge  clearly  as  the  traces  are  progressively 
summed.  A  synthetic  for  the  layer  over  a  half  space  model  is  shown  at  the  bottom. 
Again  the  correspondence  between  data  and  synthetics  suggests  that  the  composition 
of  Pg  in  terms  of  its  energy  is  much  like  that  in  the  simple  crustal  model. 
This,  however,  raises  many  interesting  questions.  The  crust  in  the  western  US. 
is  certainly  not  a  simple  layer  over  a  half  space.  It  contains  strong  gradients, 
at  least  at  the  surface  and  strong  lateral  variations  in  thickness.  This  should 
cause  large  discrepancies  between  data  and  synthetics  at  least  at  short  periods. 
We  report  our  investigations  into  these  questions  next. 

Finally,  we  presented  in  our  second  report  some  preliminary  results  indicating 
that  crustal  resonance  phases  might  be  useful  for  event  discrimination.  It 
appears  that  the  first  resonance  phase  is  primarily  a  channel  for  compressional 
energy.  The  second  is  a  much  more  efficient  channel  for  converted  energy.  At 
short  periods  in  the  real  earth,  energy  traveling  in  these  channels  will  be 
highly  scattered,  but  there  should  be  more  shear  energy  in  the  second  channel. 
Thus,  the  second  channel  may  be  more  strongly  excited  by  sources  richer  in  shear 
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PROGRESSIVE  SUM  Of  RECORDS 
vuCCS  TO  alC  OR  AMMO 


Figure  3.  Enhancement  of  crustal  resonance  phases  for  Yucca  explosions  recorded  at  ALQ  or 
ANMO  by  progressive  sum  of  seismograms.  The  synthetic  at  the  bottom  merely  indicates  that  the 
phases  have  approximately  the  right  timing  and  amplitude. 
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energy.  Figure  U  compares  averaged  three  component  records  for  Pahute  explosions 
observed  at  ANMO  to  vertical  and  radial  synthetics  for  the  simple  crustal  model. 
It  is  clear  that  the  particle  motion  changes  in  the  synthetics  between  the  first 
and  second  crustal  resonance  phase.  Because  of  interference  effects  the  second 
resonance  phase  is  suppressed  on  the  radial.  The  same  happens  in  the  stacked 
observations.  The  tangential  record  suggests  that  there  are  bursts  of  scattered 
energy  perhaps  associated  with  the  crustal  resonance  channels. 


9 


Figure  4.  Average  three-component  seismograms  recorded  at  AMMO  and  obtained  by  stacking 
the  explosions  From  Pahute  Mesa.  Evolution  of  the  particle  motion  between  the  first  and  second 
crustal  resonance  phases  is  clear.  Synthetics  were  computed  to  predict  the  observed  particle-motion 
evolution. 


FINAL  RESULTS  Ofi  SHORT  PERIOD 

Long  period  PnI  from  explosions  and  earthquakes  has  been  modeled  previously 
in  a  number  of  significant  and  remarkably  successful  investigations.  In  the 
first  of  these,  Helmberger  (1972,1973)  modeled  the  long  period  P  waveforms  from 
NTS  explosions  as  they  progressed  from  waves  trapped  in  the  crust  into  diffracted 
waves  in  the  shadow  zone  and  finally  into  waves  turning  in  the  upper  mantle. 
This  work  resulted  in  a  widely  accepted  model  for  the  velocity  and  Q  structure 
of  the  low  velocity  zone.  The  next  investigation  was  the  one  discussed  above 
by  Helmberger  and  Engen  (1980)  who  focused  on  the  waves  trapped  in  the  crustal 
wave  guide.  They  detailed  procedures  for  computing  synthetics  and  introduced 
the  term  Pnl  waves  to  describe  the  composite  of  long  period  Pn  and  PL  energy. 
In  subsequent  investigations,  Wallace  etal.  (1981)  and  Wallace  (1986)  demonstrated 
that  long  period  Pnl  waveforms  are  stable  and  deterministic  enough  to  be  utilized 
in  automated  inversion  schemes  for  source  mechanisms  or  crustal  structure.  The 
Investigations  of  Somerville  (1986),  Lefevre  and  Helmberger  (1989)  and  Holt  and 
Wallace  (1989)  demonstrated  that  the  Pnl  analysis  procedures  could  be  applied 
without  substantial  modification  to  observations  from  virtually  any  tectonic 
regime.  Figure  5  shows  the  locations  and  mechanisms  of  the  events  In  the  North 
American  continent  for  which  Pnl  waves  have  been  successfully  modeled.  The 
applicability  of  the  technique  clearly  does  not  depend  on  fault  mechanism  or  on 
tectonic  environment.  In  this  investigation,  we  have  returned  once  more  to 
observations  of  Pnl  from  NTS  explosions  as  in  the  first  studies.  However,  we 
extend  the  analysis  down  to  much  shorter  periods.  The  earlier  investigations 
all  focused  on  energy  with  periods  longer  than  about  10  seconds.  We  here  discuss 
the  analysis  of  data  with  periods  as  short  as  2  seconds. 
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Figure  5.  Location  and  focal  mechanism  of  individual  earthquakes  of  the  North  American  continent 
for  which  the  long-period  Pnl  seismograms  were  successfully  modeled. 


Additional  motivation  for  attempting  such  a  study  is  illustrated  in  Figure 
6.  There  we  compare  observed  to  synthetic  Pnl  waves  for  the  8/16/66  Caliente 
earthquake  and  the  1/9/82  New  Brunswick  earthquake.  Their  locations  are  shown 
in  Figure  5.  The  correspondence  of  the  long  period  character  in  the  data  and 
synthetics  is  very  apparent.  It  is  a  remarkable  fact  that  exactly  the  same 
crustal  structure  was  used  in  calculation  of  the  synthetics  for  these  two  events. 
It  was  the  simple  one  layer  over  a  half  space  model  introduced  in  the  Helmberger 
and  Engen  (1980)  paper.  In  fact,  this  same  simple  model  has  been  used  in  most 
of  the  above  mentioned  studies  to  model  Pnl  all  over  the  world.  The  success  of 
the  model  in  predicting  the  observations  in  the  various  tectonic  regimes  is 
somewhat  surprising  since  it  is  widely  recognized  that  the  crust  in  the  Canadian 
shield  is  almost  twice  as  thick  as  in  the  western  U.  S.  and  significantly  higher 
in  velocity.  The  general  success  of  this  simple  model  is  explained,  however, 
by  the  fact  that  at  long  period  the  earth  does  behave  like  a  layer  over  a  half 
space.  The  detailed  nature  of  the  crust  is  not  significant.  In  Figure  6,  we 
have  connected  with  short  line  segments  some  of  the  shorter  period  details  in 
the  data  and  synthetics  which  appear  to  correspond.  This  shorter  period  energy 
will  certainly  be  very  sensitive  to  details  of  crustal  structure  in  the  two 
tectonic  regions.  Yet  the  features  connected  by  the  line  segments  appear  as  if 
they  would  match  with  only  moderate  adjustments  of  the  crustal  structure. 

There  have  been  significant  advances  in  the  last  few  years  both  in  the  field 
of  seismic  data  recording  and  in  computational  seismology.  Data  is  now  recorded 
digitally  and  advanced  computers  permit  rapid  calculation  of  relatively  exact 
synthetic  seismograms  in  many  layered  structures.  In  this  study,  we  have  taken 
advantage  of  these  developments  to  try  to  model  some  of  the  finer  details  of  the 
Pnl  signals  from  NTS  explosions. 
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Figure  6.  Pnl  seismograms  from  earthquake  sources  showing  the  correlation  between  data  (dark 
lines)  and  synthetics.  Near  vertical  lines  show  the  features  that  are  well  correlated. 
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DATA  AMP  DATA  PROCESSING  Modern  broad  band  digital  recording  systems  are  rapidly 
becoming  commonplace  throughout  the  world,  particularly  near  such  important  data 
sources  as  the  Nevada  Test  Site.  Unfortunately,  methods  of  cataloguing  and 
distributing  data  remain  diverse  in  nature  and  efficiency.  Figure  7  shows  the 
eight  regional  stations  around  NTS  from  which  we  assembled  the  data  base  for 
this  investigation.  The  two  stations  ALQ  and  JAS  are  part  of  the  DWWSSN,  the 
four  stations  MNV,  ELK,  KNB  and  LAC  are  Lawrence  Livermore  National  Laboratory 
stations  and  PAS  and  PFO  are  operated  independently  by  universities.  The  data 
from  the  DWWSSN  is  the  most  generally  accessible  from  the  network.  Pnl  signals 
from  approximately  30  events  for  each  of  the  two  stations  were  plotted  and 
reviewed  for  possible  use.  Data  from  the  LLNL  network  is  not  routinely  made 
generally  available.  However,  representative  data  bases  are  occasionally 
released.  A  suite  of  2U  events  were  reviewed  for  possible  use.  The  new  broad 
band  digital  stations  at  Pinyon  Flat  and  Pasadena  were  installed  most  recently. 
The  instruments  are  new  state-of-the-art  Streckeisens ,  but  they  have  only  been 
in  operation  for  a  brief  time.  Seismic  data  from  6  nuclear  tests  were  provided 
by  the  observatory  for  possible  use  in  the  study.  The  data  from  Pasadena  station 
may  be  accessed  through  a  telephone  link,  but  only  a  brief  historical  backlog 
exists.  Signals  from  eight  recent  tests  were  collected  and  reviewed. 

The  instruments  at  the  eight  stations  vary  in  their  characteristics.  The 
DWWSSN  signals  are  comparatively  band  limited  with  respect  to  all  of  the  others. 
The  other  six  stations  record  broad  band  high  dynamic  range  signals  .  To  standardize 
our  data  base,  we  corrected  all  signals  to  the  effective  WWSSN  response.  This 
response  was  selected  because  accurate  broad  band  information  is  difficult  to 
recover  from  the  DWWSSN  data,  and  WWSSN  records  have  been  used  in  virtually  all 
of  the  previous  Pnl  studies.  Initial  review  of  the  data  indicated  that  there 
was  a  great  deal  of  stability  in  the  Pn  signals  from  the  large  events  and  also 
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consistency  in  the  PL  especially  near  the  PL  onset.  At  later  times  in  the  record, 
the  effects  of  scattering  begin  to  dominate  and  the  signals  at  any  given  station 
become  highly  variable  from  event  to  event.  Accordingly,  sections  of  the  signals 
from  all  of  the  events  being  reviewed  were  windowed  out  including  a  10  second 
leader,  all  of  the  Pn  train  and  about  15  seconds  at  the  front  of  PL.  We  will 
discuss  the  composition  of  this  last  part  of  the  signal  in  much  more  detail  in 
the  following.  We  will  show  there  that,  depending  on  choice  of  nomenclature, 
this  last  part  could  as  correctly  be  described  as  the  onset  of  the  regional  phase 
Pg  or  the  15  seconds  after  the  arrival  of  the  Moho  reflection  PmP. 

The  original  instrument  response  for  each  signal  was  deconvolved  out  and 
the  WWSSN  response  convolved  in.  A  further  bandpass  filter  was  applied  to  the 
signals.  The  form  of  the  filter  for  both  high  and  low  pass  was  a  standard  third 
order  Butterworth  given  in  the  Laplace  domain  by 

to,  oof 

/r(s)-  - i_ — - - - -  (1) 

s  ♦  tot  sz  +  cos(n/6)cots  ♦  oof 

where  s-iw.  The  above  filter  was  designed  in  the  Z-domain  using  the  bilinear 
Z- transform  and  filtering  was  performed  in  the  time  domain.  The  low  pass  value 
for  (^l,  the  cutoff  frequency,  was  set  at  0 . 6  hz .  As  discussed  in  the  next 
section,  this  made  computation  of  synthetics  using  a  wave  number  integration 
approach  relatively  practical.  The  same  Z-domain  functional  form  of  a  low  pass 
filter  becomes  a  high  pass  filter  if  s  is  replaced  by  its  inverse  and  the  cut-off 
frequency  is  defined  as  The  high  pass  cutoff  was  set  at  0.1  hr 
to  eliminate  long  period  offsets  and  drifts. 

Because  of  extant  yield  threshold  agreements,  all  events  detonated  at  NTS 
since  installation  of  the  digital  network  have  been  relatively  small  (n^  <  6) . 
It  is  well  known  that  microseismic  noise  is  a  significant  source  of  signal 
contamination  for  events  of  this  size  at  the  distances  and  period  range  under 
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consideration .  In  reviewing  the  data  base,  we  found  that  the  Pn  wavetrain  in 
particular  was  below  the  noise  level  in  very  many  cases.  We  were  able  to  find 
within  the  data  base,  however,  several  signals  for  large  events  on  quiet  days 
in  which  Pn  signal  to  noise  ratios  were  good.  Comparisons  of  these  signals 
demonstrated  that  the  waveshapes  were  consistent  through  the  Pn  and  the  onset 
of  high  amplitude  reflected  energy.  Though  the  uncertainties  associated  with 
background  noise  were  substantial,  it  was  found  that  the  entire  data  base  indicated 
that  de termini Stic  modeling  of  this  part  of  the  wavetrain  would  be  meaningful. 
Figure  8a  illustrates  this  stability  for  3  different  events  at  each  of  the  4 
LLNL  stations  for  Yucca  Flat  events.  The  move  out  of  Pn  in  front  of  the  strong 
reflected  energy  as  distance  increases  is  apparent.  The  consistency  of  the  short 
period  details  in  the  waveforms  is  even  more  striking  than  in  the  earthquake 
data  in  Figure  6.  A  similar  comparison  of  signals  for  Pahute  Mesa  events  is 
presented  in  Figure  8b.  The  strong  correlation  of  the  signals  is  again  apparent. 
There  are  differences  between  the  Yucca  Flat  and  Pahute  Mesa  signals  at  the  same 
station,  but  this  is  not  surprising  given  the  variation  in  distance  to  the  sites. 

Figure  9a  compares  the  waveforms  for  3  of  the  remaining  stations  in  the 
digital  net  for  Yucca  Flat  events.  Because  of  the  very  limited  number  of  events 
available  from  PAS,  we  were  only  able  to  find  two  events  with  acceptable  signal 
to  noise.  No  suitable  two  Yucca  Flat  events  at  all  were  in  the  6  signals  available 
from  PFO.  Nonetheless,  our  basic  observation  regarding  the  stability  of  the 
beginning  of  the  wave  train  holds  true.  The  observations  from  ALQ  are  particularly 
remarkable  given  that  the  station  is  so  far  away.  The  Pn  wave  train  is  almost 
30  seconds  long,  but  the  strong  reflected  energy  remains  very  stable.  A 
deterministic  modeling  study  would  clearly  be  appropriate.  Figure  9b  presents 


18 


Figure  8.  Recorded  Pnl  seismograms  from  both  Yucca  flat  and  Pahute  Mesa  on  the  LLNL  network 
after  convolving  with  a  WWSSN  LP  instrument  and  a  recursive  bandpass  filter  (low-cut-  0.1  Hz. 
high  cut  -  0.6  Hz  and  third-order  Butterworih  causal),  (a)  For  Yucca  flat  explosions  and  (b)  for 
Pahute  Mesa  explosions.  Data  show  remarkable  stability.  Original  instrument  response  was 
deconvolved. 
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a  similar  comparison  for  Pahute  Mesa  events  for  the  remaining  4  stations  in  the 
digital  net  including  PFO  this  time.  Modeling  these  waveforms  has  been  the 
primary  goal  of  this  study. 

Numerical  Methods  Pnl  waveforms  can  be  synthesized  at  regional  distances  by  a 
variety  of  techniques.  The  two  most  practical  are  generalized  rays  and  full-wave 
synthesis.  The  generalized  ray  approach  is  the  most  efficient  for  a  layer  over 
a  half  space  model  and  has  been  used  in  most  previous  Pni  studies.  It  becomes 
difficult  to  obtain  exact  results  in  a  multilayered  model  because  the  ray  sum 
becomes  intractably  large.  One  real  utility  of  the  method  is  that  it  allows 
decomposition  of  the  wave  field  into  individual  rays.  We  will  take  advantage 
of  this  facet  to  some  extent  in  the  following.  The  full-wave  synthesis  approach 
can  be  effectively  used  to  produce  an  accurate  estimate  of  the  wave  field  in  a 
finite  frequency  band  for  a  complex  layered  mode1  However,  some  steps  must  be 
taken  to  band  limit  both  the  waven”niber  and  frequency  integrals  to  keep  the 
method  truly  practical.  Alternate  approaches  to  computing  the  synthetic  seis¬ 
mograms  are  the  frequency-wavenumber  (f-k)  scheme  (Wang  and  Herrmann,  1980; 
Bouchon,  1981)  and  the  reflectivity  approach  (Kennett,  1980;  Mallick  and  Frazer, 
1987).  Our  primary  numerical  tool  in  this  work  has  been  a  modified  version  of 
the  Wang  and  Herrmann  (1980)  method. 

We  developed  a  vectorized  computer  program  based  on  the  formulation  presented 
by  Wang  and  Herrmann  (1980).  For  a  medium  consisting  of  n  laterally  homogeneous 
layers  and  a  source  embedded  at  an  interface  a  the  f-k  kernels  are  evaluated 
using  the  expression 


I'A  l  (  K.x  \“Z„  \ 
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(the  a,/s  being  the  elements  of  the  original  Haskell's  (1964)  matrix).  The 

explicit  forms  of  F E«,  a  and  (  the  source  coefficients)  are  given  in  Wang 
and  Herrmann  (1980) .  Expression  (2)  together  with  the  expressions  (3)  allow  the 
propagation  of  the  stress-displacement  vector  from  a  half-space  through  a  layered 
medium  with  a  source  discontinuity  to  the  free  surface.  It  involves  propagation 
of  the  stress-displacement  vector  with  the  use  of  compound  matrices  R  |  jf  and  .V  I \f 
(Dunkin,  1965;  Watson,  1970).  These  compound  matrices  are  the  product  of  matrices 
of  subdeterminants.  They  have  an  advantage  over  the  Haskell  matrices  in  that 
the  square  of  the  exponential  terms  have  been  canceled  out  in  the  analytical 
formulation  of  a  |^(  matrix.  This  eliminates  the  loss  of  significant  figures  in 
the  multiplication  and  subtraction  of  the  original  Haskell  matrix  elements 
(Dunkin,  1965) . 

The  f-k  kernels  U(r .  k. ou) are  integrated  to  obtain  the  medium  response  L'{r.uo) 
using  the  Fourier-Bessel  transformation 

I/(r.u>)-  J  L’(.r.k.M)kJ„{kr)dk  (4) 

0 

For  the  evaluation  of  the  regional  wavefield,  we  have  approximated  the  Bessel 
functions  by  the  principal  asymptotic  form  of  the  Hankel  function  representation 
of  the  out-going  wavefield.  With  these  approximations,  the  final  expression 
used  to  integrate  the  kernels  is 
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We  use  f(uj)*  /  f(t)e~"“'dtfoT  FFT.  Wavenumber  Integrals  are  evaluated  using  equal 


sampling  in  k  domain  (Bouchon,  1981)  and  inversion  to  the  time  domain  using  an 
FFT  algorithm.  The  code  and  the  reflectivity  code  developed  by  Mallick  and 
Frazer  (1988),  based  on  Kennett  (1980),  were  calibrated  against  each  other.  The 
results  agreed  very  closely,  and  the  f-k  algorithm  proved  to  be  more  economical 
for  computing  the  free  surface  displacement  in  the  frequency  band  of  our  interest. 

The  ' — 'yw  of  interest  in  Pnl  studies  is  approximately  one  minute  after 
the  first  arrival.  At  regional  ranges,  the  entire  signal  is  dispersed  considerably 
lasting  for  up  to  10  minutes.  Since  the  spectrum  of  the  complete  signal  is 
automatically  computed  by  the  method,  this  later  energy  can  strongly  affect  the 
spectrum  in  the  time  window  of  interest  through  the  phenomenon  of  wrap  around 
if  a  short  time  window  is  chosen.  Several  possible  approaches  for  overcoming 
this  problem  have  been  suggested  in  the  literature,  and  we  have  investigated 
three  of  them.  These  include  integration  on  a  displaced  contour  in  the  frequency 
domain  (Phinney,  1965;  Bouchon  and  Aki ,  1977;  Spudich  and  Ascher ,  1983),  wavenumber 
filtering  and  computation  of  an  extended  time  series  to  damp  out  wrap  around, 
and  we  have  tested  all  three.  We  have  also  investigated  the  importance  of  Q. 

The  damping  technique  of  Phinney  (1965)  involves  giving  frequency  a  small, 
constant  imaginary  part  a  which  has  a  two- fold  purpose.  One  is  the  shifting  of 
the  poles  of  the  f-k  kernels  from  the  real  k  axis  into  the  fourth  quadrant  of 
the  complex  k  plane.  Consequently,  integration  can  be  performed  along  the  real 
k  axis  without  any  problem.  The  use  of  complex  frequency  also  has  the  effect 
of  multiplying  the  true  response  f(t)  by  a  factor  of  by  e'at  Mallick  and  Frazer 
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(1987)  showed  via  an  elegant  Fourier  transform  pair  that  the  evaluation  of 
response  f(c)  using  the  complex  frequency  technique  produces  a  time  domain 
response  equivalent  to 


e -*'[  fit)*  Z  fit*nT)e 
L  «*o 


-2.-  1 ,  1 .2.... 


(6) 


i.e.  it  suppresses  the  response  £(t+nT)  outside  the  time  window  T  by  a  factor 
e-a<t+nT)  They  suggest  a  rule  of  thumb  in  selecting  a  of  l/D-e'aT  where  the 
maximum  allowable  value  of  D  is  50  and  T  is  the  total  length  of  the  time  series. 
This  damping  is  then  corrected  for  in  the  final  time  domain  synthetics.  Bouchon 
(1979)  suggested  that  a  should  be  chosen  such  that  e‘ar  «  1.  Herrmann  (1985) 
pointed  out  that  aT-2  is  a  reasonable  choice.  He  also  suggested  that  some  trial 
and  error  testing  may  be  necessary  to  set  this  parameter.  As  we  show  in  the 
following,  our  experiments  with  a  simple  crustal  model  indicate  that  a  higher 
value  of  D  will  cause  spurious  noise  which  originates  due  to  the  values  of  u 
which  no  longer  satisfy  the  condition  w»o. 

Our  test  model  was  simply  a  layer  over  a  half  space  with  a  few  low  velocity 
layers  on  top  to  represent  the  surface  gradient.  Our  initial  selection  for  time 
series  length  was  1024  points  with  a  Nyquist  frequency  of  5  hz.  D  was  set  at 
50.  The  source  is  a  shallow  explosion  and  the  range  is  300  km.  Here  and 
throughout  we  utilize  the  Murphy  (1971)  wet  tuff  source  time  history.  The  results 
are  shown  in  the  top  row  of  Figure  10.  The  notation  above  the  trace  indicates 
the  value  of  a,  that  no  Q  correction  is  applied  and  no  wave  number  filter  is 
used.  The  notation  is  self  explanatory  for  the  rest  of  the  traces.  Though  the 
synthetic  looks  generally  reasonable  the  contamination  of  the  Pnl  is  severe.  The 
short  trace  on  the  right  is  an  amplification  of  the  Pnl  from  the  left  trace,  and 
it  clearly  shows  high  amplitude  acausal  energy  at  the  beginning.  We  next  attempted 
increasing  the  number  of  points  to  4096.  We  reduced  the  Nyquist  to  2.5  in  an 
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Figure  10.  A  suite  of  theoretical  WWSSN  long-period  seismograms  computed  using  a  vectorized 
f-k  algorithm  from  a  half-space  model  with  a  gradient-like  structure  in  the  upper  six  kilometers. 
Seismograms  on  left  demonstrate  the  effects  of  various  Q,  ray-parameter  filters,  complex  frequencies 
and  length  of  time  series.  On  right  are  plotted  the  Pn\  seismograms  of  this  study. 


26 


effort  to  Improve  efficiency  and  reduced  D  to  10  to  improve  accuracy.  The 
required  increase  in  computation  time  for  the  longer  time  series  was  still  very 
large.  The  traces  in  the  second  row  in  Figure  10  show  the  result.  The  complete 
signal  on  the  left  shows  unrealistically  high  amplitude  higher  mode  surface 
waves.  In  the  real  earth  this  type  of  energy  would  ordinarily  be  attenuated  or 
scattered  out.  Anelasticity  was  introduced  through  complex  compressional  and 
shear  wave  velocities  (a,. 0,)  complex  as  given  by 


a,(u>) 


P,(u>) 


(7) 


where  wr  is  reference  frequency  (Aki  and  Richards,  1980).  We  chose  this  reference 
frequency  to  be  1  Hz.  A  relationship  of  Qa-2Q^  was  used  for  Qa.  The  near- surface 
Qp  structure  was  taken  from  the  basin  and  range  province  of  the  Western  United 
States  (Wen,  1989).  The  realistic  looking  results  are  shown  in  the  third  row. 
The  Pnl  portions  of  the  seismogram  on  the  right  of  rows  two  and  three  can  be 
considered  as  the  "correct"  synthetics.  Our  additional  efforts  were  directed 
at  finding  a  way  to  compute  them  efficiently.  A  comparison  between  the  Pnl  waves 
in  rows  2  and  3  indicates  that  attenuation  has  little  effect  on  the  waveshape 
and  a  small  effect  on  amplitude.  Thus,  it  is  not  important  to  know  the  Q  structure 
accurately  to  compute  Pnl  in  this  frequency  band.  We  next  considered  the 
possibility  of  wave  number  filtering  and  used  the  same  Q  structure. 

The  filter  we  applied  was  a  simple  cosine  taper  defined  in  ray  parameter 
domain  by 


2^  1  _cos 
F(P)~{  10 
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The  values  of  p^,,  and  px  were  fixed  at  0.0001  and  0.033  sec  km'1  ,  and  Che  effects 
of  the  high  p  cutoff  were  tested.  Figure  10  shows  the  results  of  p2-0-285  and 
pmmx-0.288  sec  km'1  (WF1),  p2-0.222  and  P...-0 . 244  sec  km'1  (WF2)  and  p2-0.200  and 
P-..-0  ■  220  sec  km'1  (WF3)  .  The  goal  is  to  select  parameters  which  preserve  both 
the  shape  and  amplitude  of  the  ?nl  waveform.  As  indicated  in  the  figure,  WF1 
leaves  too  much  long  period  and  WF3  distorts  the  waveform  substantially.  The 
parameters  of  WF2  give  acceptable  results.  We  examined  the  possibility  of 
abandoning  the  complex  frequency  damping  in  light  of  the  efficiency  of  the 
wavenumber  filter  in  reducing  wrap  around.  As  shown  in  the  bottom  row,  some 
damping  is  required.  To  summarize ,  we  found  that  for  the  types  of  models  discussed 
in  the  following  the  optimal  way  to  compute  Pnl  synthetics  is  to  use  a  time  series 
of  1024  points,  a  Nyquist  frequency  of  2.5  hz ,  WF2  and  an  a  value  of  .0028  hz . 
Numerical  Experiments  Figures  4  and  5  demonstrate  that  the  short  period  waveshape 
of  Pnl  is  relatively  stable  up  through  the  first  ten  seconds  after  the  onset  of 
high  amplitude  signal  after  P„.  We  discuss  here  a  few  numerical  experiments 
designed  to  elucidate  the  composition  of  this  signal  and  its  sensitivity  to 
details  of  the  crustal  structure.  We  begin  with  generalized  ray  calculations 
in  the  Helmberger  and  Engen  (1980)  layer  over  a  half  space  crustal  model.  We 
have  computed  a  section  of  displacement  Green's  functions  for  an  isotropic 
compressional  source  and  show  the  result  on  the  left  of  Figure  11.  We  included 
254  rays  in  the  calculation  corresponding  to  all  paths  with  three  or  less 
reverberations  in  the  crustal  layer.  Helmberger  and  Engen  (1980)  show  the  build 
up  of  the  responses  as  up  to  4094  generalized  rays  are  successively  added  in 
corresponding  to  5  reverberations  in  the  crust,  and  generally  speaking  the  next 
two  levels  of  reverberation  do  not  contribute  significantly  to  the  time  window 


of  interest  here . 
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Figure  11.  Broadband  seismograms  Co  show  amplitude  responses  of  individual 
rays  compuced  by  GENRAY88  (generalized  ray  theory)  using  a  layer  over  half- 
space  model  of  Helmberger  and  Engen  (1980).  Major  contributing  generalized 
rays  are  identified  in  the  suite  of  right  most  seismograms. 


29 


p 


The  high  amplitude  portion  of  the  responses  appears  to  consist  of  a  rhythmic 
sequence  of  sharp  pulses  which,  judging  from  their  complexity  of  shape,  are 
probably  the  result  of  several  different  generalized  rays.  We  have  examined  the 
response  to  small  subsets  of  rays  to  attempt  to  identify  which  paths  contribute 
the  primary  energy  to  these  pulses.  We  found,  as  might  have  been  expected,  that 
the  rhythmic  sequence  is  related  to  successive  orders  of  crustal  reverberation 
with  one  relatively  interesting  complexity.  A  response  section  built  from  only 
9  rays  which  appear  to  contribute  most  to  the  energy  to  the  beginning  of  the 
Green's  function  is  shown  on  the  right  of  the  figure.  As  indicated  there,  the 
first  high  amplitude  pulse  is  made  up  of  the  post  critical  mantle  reflection 
PmP,  the  direct  ray  noted  as  Pd  and  pPraP.  The  step- like  Pn  wave  is  identified 
separately,  but  in  a  generalized  ray  sense  it  is  part  of  the  PmP  ray.  The  second 
reverberation  of  PmP  is  indicated  in  the  diagram  as  2PmP.  It  does  not  contribute 
significantly  to  the  second  high  amplitude  resonance  pulse.  This  pulse  is 
composed  primarily  by  PmPSmP  and  PmSPmP.  In  other  words,  the  second  high  amplitude 
sp-ke  is  caused  by  rays  with  at  least  one  converted  leg.  The  pSmS  and  PmPSmS 
follow  this  high  amplitude  spike.  The  PmS  phase  precedes  the  spike.  The  later 
resonances  seen  in  the  complete  responses  on  the  left  of  the  figure  are  due  to 
similar  complex  groups  cf  rays  which  arrive  closely  in  time.  Some  of  the  results 
regarding  the  timing  and  amplitude  associated  with  the  contributions  of  certain 
converted  rays  in  a  layer  over  a  half  space  are  counter  intuitive  because  the 
rays  follow  very  asymmetric  paths.  S  legs  are  almost  vertical  and  P  legs  are 
relatively  flat.  The  amplitudes  of  the  rays  are  controlled  by  the  ranges  at 
which  critical  angle  occurs. 

There  are  two  points  we  intend  to  make  through  this  exercise.  The  first 
is  that  factors  which  might  tend  to  change  conversion  coefficients  and  travel 
times  of  converted  phases  should  be  important  in  understanding  the  onset  of  high 
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amplitude  signal  we  are  crying  to  model  in  this  work.  The  second  is  that  the 
very  beginning  of  this  signal  may  properly  be  referred  to  as  PmP,  though  the 
details  of  propagation  of  this  energy  will  be  much  more  complex  in  a  realistic 
crust  model  than  in  a  layer  over  a  half  space. 

To  begin  our  investigations  into  the  behavior  of  Pnl,  we  designed  several 
simple  test  models  which  had  some  of  the  major  features  one  would  expect  to  see 
on  average  in  a  real  crust.  They  are  shown  in  Figure  12.  The  first  5  models 
are  of  a  surface  gradient.  Model  1  has  a  velocity  decrease  of  about  20%  in  P 
and  30%  in  S  over  a  range  of  6  km  in  depth  on  top  of  the  standard  one  layer 
model.  The  Green's  functions  for  the  one  layer  model  and  Model  1  are  compared 
in  the  left  two  columns  of  Figure  13.  The  ranges  shown  are  representative  of 
those  at  which  we  have  actual  observations.  The  effect  of  the  simple  gradient 
is  very  pronounced  at  high  frequency.  The  response  with  the  surface  gradient 
present  is  much  more  complex  and  apparently  higher  in  frequency.  The  primary 
reasons  are  that  the  free  surface  reflection  coefficients  change  substantially 
along  with  the  timing  of  rays.  In  the  one  layer  model,  the  conversion  of  P  to 
S  energy  (rays  like  pS)  is  much  stronger,  whereas  in  Model  1  more  energy  stays 
in  the  compressional  mode  (rays  like  pP) .  In  model  2,  we  experimented  with  an 
additional  substantial  drop  in  the  P  velocity  with  the  S  velocity  fixed.  In 
model  3,  we  did  the  reverse.  The  Green's  functions  did  not  change  substantially 
in  shape  as  illustrated  in  the  right  two  columns  of  Figure  13.  We  therefore 
infer  that  the  waveform  results  are  not  strongly  sensitive  to  the  ratio  of  P  to 
S  velocity  in  the  crust  gradient.  Amplitudes  vary  by  about  a  factor  of  2.  Models 
and  5  represent  progressively  stronger  gradients  in  both  P  and  S  and  will  be 
discussed  further  in  the  following.  The  Green’s  functions  shown  in  Figure  13 
were  computed  over  a  frequency  range  extending  up  to  5  hz .  The  data  we  wish  to 
model  has  been  filtered  and  convolved  with  the  long  period  WWSSN  instrument 
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Figure  12.  Canonical  crustal  models.  Top  two  panels  show  the  variation  of  P  and  S  wave  e  o 
within  the  upper  8  km  of  the  crust  (models  extends  up  to  32  Km.  Half-space  P  and  s  ^ocit 
are  8.2  km  sec  and  4.5  km/sec).  The  lower  two  panels  show  the  structure  across  the 
discontinuity.  Models  1-3  are  for  the  perturbation  study  discussed  in  the  text. 
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CRUST  MODEL  RESPONSES 


Figure  13.  Broadband  f-k  seismograms  with  a  Nyquist  of  5  Hz  computed  for  models  1-3  including 
a  layer  over  half-space  model.  All  seismograms  were  convolved  with  a  Mueller-Murphy  source. 
The  effects  of  the  upper  crust  can  be  clearly  noticed  in  the  change  of  the  frequencies  and  complexity 
of  the  waveforms. 
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response.  It  contains  little  energy  at  frequencies  higher  than  1  hz .  To  confirm 
that  it  is  important  to  account  for  the  surface  gradient  in  our  modeling,  we 
convolved  the  Green's  functions  with  the  instrument  response.  The  results  are 
shown  in  Figure  14.  The  arrangement  of  the  columns  is  the  same  as  in  Figure  13. 
The  surface  gradient  continues  to  be  important  in  the  lower  frequency  band. 

We  next  examined  the  effects  of  a  gradient  at  the  crust-mantle  interface. 
We  computed  synthetics  at  the  same  ranges  for  models  6  and  7  in  Figure  12.  The 
velocities  above  and  below  the  interface  are  assumed  to  be  constant.  We  found 
that  the  effects  of  the  gradients  were  not  significant  at  the  level  of  the  scatter 
in  the  waveform  data  and  certainly  not  nearly  as  important  as  the  surface  gradient. 
The  Green's  functions  for  models  6  and  7  are  compared  to  those  for  the  one  layer 
model  in  Figure  15  and  with  the  WWSSN  response  in  Figure  16.  A  gradient  at  the 
Moho  must  certainly  have  significant  effects  on  waveforms  at  some  ranges, 
particularly  when  the  first  arrival  is  triplicating  through  it.  This  occurs  at 
around  150  km  in  the  western  U.S.  The  closest  range  from  which  we  have  data  is 
237  km,  so  all  of  the  energy  near  the  onset  of  Pg  is  well  past  critical  angle 
and  strongly  trapped  in  the  crust.  This  explains  its  insensitivity  to  the  details 
of  the  crust-mantle  transition  zone.  The  results  of  these  numerical  experiments 
also  depend  on  how  realistic  the  monels  we  tested  are.  In  fact,  we  selected 
them  based  on  previous  estimates  of  the  strength  of  the  surface  and  Moho  gradients . 
Figure  17  compares  our  model  4  to  the  Pahute  Mesa  crust  model  of  Hartzell  et  al. 
(1983).  In  fact  they  are  very  similar,  and,  as  we  shall  show  in  the  following, 
this  model  is  very  consistent  with  the  data.  On  the  right  of  the  figure,  we 
show  a  model  developed  for  the  crust-mantle  transition  for  the  path  from  NTS  to 
Las  Vegas.  This  work  is  discussed  in  the  final  section  of  this  report.  This 
realistic  gradient  model  is  very  similar  to  our  model  7. 
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Figure  14.  Same  seismograms  as  shown  in  Figure  13,  but  after  convolving  with  a  long-period 
WWSSN  instrument. 
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Figure  15.  Broadband  f-k  seismograms  with  a  Mueller-Murphy  source  function  computed  using 
a  transition  zone  3  km  thick  crust  across  the  Moho  discontinuity  in  the  Helmberger  and  Engen 
(1980)  crustal  model.  The  effect  can  be  seen  in  the  changes  of  the  relative  amplitude  of  phases. 
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Figure  16.  Same  seismograms  as  shown  in  Figure  16,  but  after  convolving  with  a  long-period 
WWSSN  instrument.  The  effect  of  the  transition  zone  is  small. 
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Figure  17.  Crustal  Models  with  both  P  and  S  wave  velocities.  On  the  left  are  our  preferred  model 
4  and  Hartzell  et  at  (1982,  dark  lines)  and  on  the  right  is  a  crustal  model  developed  from  NTS  to 
Las  Vegas. 
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MODELING  RESULTS  We  proceeded  then  with  the  assumption  that  the  most  important 
portion  of  the  crust  model  to  adjust  was  the  lid  gradient  in  order  to  improve 
the  fit  to  short  period  Pnl.  Since  our  numerical  experiments  indicated  that  the 
ratio  of  P  to  S  velocities  was  not  important,  we  tested  the  effects  of  increasingly 
intense  crustal  gradients.  We  refer  to  the  layer  over  a  half  space  model  as 
model  0  and  the  increasing  gradients  as  models  1,  4  and  5  as  illustrated  in 
Figure  12.  The  results  of  the  experiments  are  shown  in  Figures  18  through  21. 
For  each  of  the  eight  stations,  the  data  is  shown  at  the  top  and  the  synthetics 
for  each  of  the  four  models  below. 

As  our  numerical  experiments  would  have  indicated  (see  Figure  13) ,  addition 
of  the  surface  gradient  causes  a  notable  increase  in  the  amount  of  high  frequency 
energy.  Model  0,  the  layer  over  a  half  space,  is  obviously  too  long  in  period 
to  satisfy  the  observations.  This  is  particularly  apparent  in  the  first  arriving 
energy  or  Pn.  The  effect  of  the  weakest  gradient  model,  model  1,  is  not  strong 
enough  and  the  synthetics  are  still  too  low  frequency.  The  gradient  model  4 
predicts  approximately  the  correct  frequency  content,  and  its  synthetics  are 
shown  directly  beneath  the  data.  It  provides  some  remarkable  waveform  fits  to 
Pn  and  the  onset  of  PL.  The  model  4  results  are  shown  by  themselves  in  Figures 
22  and  23.  Data  is  shown  in  dark  line  over  the  synthetics  below.  Some  of  the 
more  notable  fits  are  to  ELK  from  JORNADA,  KNB  from  GIBNE,  ALQ  from  KAPPELI  and 
JAS  from  JORNADA.  It  is  important  to  emphasize  that  a  single  vertically  stratified 
model  has  been  utilized  to  fit  all  of  the  data.  No  doubt  strong  lateral  variations 
exist,  and  we  could  have  obtained  a  better  fit  by  using  a  separate  model  for 
each  path.  However,  our  main  point  here  is  to  prove  the  efficacy  of  using 
relatively  simple  models  to  predict  the  beginning  of  Pnl  down  to  periods  as  short 
as  2  seconds . 
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Figure  18.  A  comparison  between  data  (top  seismogram)  and  synthetic  seismograms  computed 
from  four  canonical  crustal  model  The  comparison  is  made  for  Yucca  flat  explosions  at  f>  ur 
LLNL  stations. 
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Figure  19.  A  comparison  between  data  (top  seismogram)  and  synthetic  seismograms  computed 
from  four  canonical  crustal  model  The  comparison  is  made  for  Pahute  Mesa  explosions  at  four 
LLNL.  stations. 
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Figure  20.  A  comparison  between  data  (top  seismogram)  and  synthetic  seismo¬ 
grams  computed  from  four  canonical  crustal  model.  Comparison  is  made  for 
Yucca  flat  explosions  at  four  stations-  JAS,  PAS,  PFO  and  ALQ. 
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Figure  21.  A  comparison  between  data  (top  seismogram)  and  synthetic  seismo¬ 
grams  computed  from  four  canonical  crustal  model.  Comparison  is  made  for 
Pahute  Mesa  explosions  at  three  stations  -  JAS ,  PAS  and  ALQ. 


Figure  22.  A  comparison  between  the  data  (top  seismograms)  both  from  Yucca  flat  and  Pahute 
Mesa  and  synthetic  seismograms  computed  from  our  preferred  model  4  The  initial  part  of  the 
synthetic  seismograms  shows  remarkable  correlation  with  data.  The  comparison  is  shown  for  the 
four  LLNL  stations. 
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Figure  23  A  comparison  between  the  data  (top  seismograms)  both  from  Yucca  tlat  and  Pahuu 
Mesa  and  synthetic  seismograms  computed  from  our  preferred  model  4  Initial  part  ofthesynthet. 
seismograms  shows  remarkable  correlation  with  data.  The  comparison  is  shown  for  the  station 

-  JAS.  PAS,  PFO  and  ALQ. 


I2UXIAL  VM  Q S  £q  WAYEfQKl  DISCRIMINANTS 

The  most  successful  regional  discriminant  developed  to  date  is  based  on 
systematic  differences  in  the  spectral  content  of  explosion  and  earthquake  sources 
(Murphy  and  Bennett,  1982;  Taylor  et  al.  ,  1988).  Perhaps  surprisingly,  but 
nonetheless  consistently,  NTS  explosions  are  reduced  in  levels  of  high  frequency 
energy  with  respect  to  earthquakes  of  the  same  magnitude.  There  is,  as  yet,  no 
physical  explanation  for  this  difference,  so  in  spite  of  the  success  of  the 
approach  in  the  western  U.S.  there  are  serious  questions  to  be  answered  regarding 
its  applicability  in  other  regions  of  the  world.  In  reports  two  and  three,  we 
discussed  the  development  of  another  somewhat  less  accurate,  but  more  physically 
based  discriminant  which  relies  on  the  waveform  of  regional  Pn  waves.  We  will 
summarize  the  key  results  here. 

Because  of  the  results  to  be  presented  in  the  following  section,  we  begin 
by  emphasizing  that  the  development  of  the  regional  waveform  discriminant  was 
rooted  in  past  experience  with  waveforms  of  teleselsmlc  P  waves  from  nuclear 
explosions.  Figure  24  is  taken  from  Burdick  et  al .  (1984).  It  shows  observed 
and  synthetic,  long  and  short  period  records  from  the  nuclear  test,  CANNIKIN. 
Arrows  on  the  left  of  the  figure  draw  attention  to  a  subtle  feature  in  the  short 
period  records  which  is  associated  with  the  arrival  of  the  phase  pP.  The  records 
where  the  feature  appears  are  assumed  to  be  along  high  Q  paths  (t*  -  0.8  s). 
Along  lower  Q  paths  the  feature  washes  out.  The  long  periods  are  not  affected 
by  pP  in  a  clearly  visible  way. 

This  previous  experience  with  the  effects  of  pP  on  teleseismic  short  period 
P  waves  proved  valuable  in  interpreting  regional  Pn  waveforms.  When  the  digital 
signals  from  NTS  explosions  were  averaged  (as  in  Figures  3  and  4)  ,  it  was  observed 
that  the  average  Pn  waveform  was  very  similar  to  those  on  the  left  of  Figure  24 
although  much  shorter  in  period.  Average  waveforms  at  5  stations  in  the  digital 
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Figure  24.  Comparison  of  observed  short  and  long  period  P  waves  with  synthetic  for  a  range  of 
t *  values  at  teleseismic  distances.  As  t •  increases,  the  interference  in  the  second  upswing  produced 
by  pP  becomes  less  apparent.  pP  can  be  observed  as  a  distinctive  double  peak  and  is  produced  in 
the  synthetics  with  an  artificial  delay. 
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net  are  shown  in  Figure  25.  The  consistent  splitting  of  the  second  upswing  is 
indicated  by  the  arrows.  The  clear  implication  is  that  the  physics  of  the  wave 
propagation  of  short  period  teleseismic  P  and  short  period  Pn  is  very  comparable. 
More  precisely,  the  interaction  of  pP  with  P  must  be  similar  in  the  two  instances. 
This  is  not  unreasonable  in  that  the  apparent  velocity  of  Pn  is  about  8  km/s  and 
that  of  teleseismic  P  only  increases  to  about  12  km/s  at  30°.  The  associated 
change  in  pP  timing  and  amplitude  is  small.  The  implications  for  regional 
discrimination  are  clear.  Only  very  shallow  sources  like  explosions  will  have 
depth  phases  at  very  short  times.  Earthquake  depth  phases  will  be  much  later. 
To  test  the  performance  of  this  discriminant,  we  assembled  a  set  of  Pn  waveforms 
from  small  earthquakes  near  NTS,  windowed  out  the  first  three  seconds  of  Pn  and 
measured  the  correlation  with  the  average  Pn  waveforms  of  explosions  like  those 
shown  in  Figure  25.  A  similar  procedure  was  carried  out  on  the  explosion  data 
base.  The  results  from  one  station  (JAS)  are  shown  in  Figure  26.  The  explosions 
are  displayed  as  stars  and  the  earthquakes  as  crosses.  The  separation  of  the 
populations  is  good  enough  to  warrant  more  study  of  this  discriminant. 

The  discrimination  capability  illustrated  in  Figure  26  only  demonstrates 
that  the  Pn  waveform  of  explosions  is  stable  at  JAS  and  consistently  different 
from  earthquake  waveforms.  However,  the  similarity  of  the  waveforms  in  Figure 
25  suggests  that  the  shape  of  the  explosion  waveform  is  consistent  from  station 
to  station.  That  this  is  indeed  the  case  is  shown  in  Figure  27  where  the  average 
JAS  explosion  waveform  has  been  correlated  with  the  explosion  data  base  at  MNV. 
Moderately  good  event  discrimination  is  still  possible.  These  studies  suggest 
that  one  method  to  achieve  discrimination  with  a  regional  net  will  be  to 
continuously  monitor  the  waveform  of  Pn  onset  at  all  of  the  stations  and  to  test 
in  an  objective  fashion  for  the  presence  of  depth  phases. 
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Figure  25.  The  average  Pn  waveform  for  Pahute  explosions  observed  at  stations  in  the  western 
U.S.  digital  net.  The  split  third  swing  is  indicated  by  arrows. 
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Figure  26.  Discrimination  of  Yucca  explosions  from  earthquakes  using  correlation  with  the  average 
Pn  waveform  at  JAS.  The  explosions  are  stars  and  the  earthquakes  are  crosses. 
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Figure  27.  Discrimination  of  Yucca  explosions  from  earthquakes  using  the  average  Pn  waveform 
from  MNV  and  the  data  base  from  JAS.  Explosions  are  stars  and  earthquake  are  crosses.  The 
separation  of  the  populations  is  comparable  to  that  in  Figure  26  illustrating  the  transportability  of 
the  discriminant. 
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One  final  study  from  report  3  which  needs  to  be  reviewed  here  is  a  forward 
modeling  study  of  the  broad  band  Pn  waveforms  from  the  LLNL  stations  which  was 
carried  out.  A  typical  result  is  shown  in  Figure  28  where  the  waveforms  of  four 
Pahute  events  as  recorded  at  MNV  are  displayed.  The  data  are  shown  as  solid 
lines  and  synthetics  for  an  appropriate  crustal  model  and  explosion  source  are 
shown  as  dashed  lines.  Because  the  instruments  are  broad  band,  deconvolving  the 
response  out  is  a  stable  operation  which  has  been  carried  out.  The  traces  shown 
are  true  ground  velocity.  The  arrows  indicate  the  arrival  of  pPn  in  the  data 
and  synthetics.  The  observed  pPn  is  consistently  much  later  than  the  elastic 
predictions.  This  discrepancy  was  observed  for  all  events  at  all  stations  where 
the  waveforms  could  be  modeled.  A  similar  result  has  been  reported  in  most 
studies  where  pP  times  from  nuclear  tests  have  been  measured.  These  results  are 
relevant  to  the  new  work  discussed  next.  More  details  of  the  investigation  are 
available  in  the  third  report  (Burdick  et  al.,  1988). 
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Eq  WAVEEORH  DISCRIMINANTS  Aim  SPALL 

The  work  in  the  final  quarter  of  this  project  on  the  Pn  waveform  discriminant 
has  been  directed  at  developing  a  further  appreciation  of  the  physics  behind  it. 
We  have  examined  the  conclusion  from  report  3  that  the  Pn  at  regional  distances 
appears  more  like  a  turning  ray  in  a  smooth  positive  gradient  than  a  true  head 
wave  on  the  crust  mantle  Interface.  To  accomplish  this,  we  have  modeled  the 
transition  of  Pn  from  a  second  to  a  first  arrival.  We  found  that  the  Moho  under 
NTS  behaves  much  more  like  a  smoothly  graded  structure  than  a  sharp  discontinuity . 
We  have  also  developed  a  physical  explanation  for  the  apparent  delay  of  pPn  as 
illustrated  in  Figure  28.  We  have  also  examined  the  possible  role  of  spall  in 
explaining  the  anomalous  late  arrival  of  pPn. 

We  have  already  reported  on  interpretation  of  the  majority  of  the  available 
digitally  recorded  data  from  NTS.  It  would,  of  course,  be  desirable  to  analyze 
even  more.  After  some  investigation,  we  learned  that  regional  data  had  been 
recorded  for  several  decades  by  DOE  on  analog  tape  and  a  portion  of  it  had  been 
digitized  with  a  high  speed  A  to  D  converter.  Figure  29  shows  the  locations  of 
the  stations  which  were  operated  for  an  extended  period  of  time.  Most  of  the 
instruments  were  L7s  which  are  very  broad  band  and,  given  the  frequency  of  NTS 
testing,  each  recorded  a  substantial  data  base.  Unfortunately,  most  of  the 
stations  are  closer  than  Pn  cross  over  distance,  so  they  can  not  be  used  to  test 
the  waveform  discriminant  directly.  They  do  provide  a  significant  data  base  for 
investigating  other  types  of  discriminants  which  are  applicable  at  shorter  ranges. 
The  stations  at  Las  Vegas  are  the  most  distant,  and  they  are  slightly  beyond 
cross  over  range.  The  top  of  Figure  30  shows  an  enlarged  map  of  the  stations, 
and  the  bottom  shows  a  seismic  section  of  the  P  waves  observed  from  the  Yucca 
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Figure  28.  Observed  versus  synthetic  Pn  waveforms  at  MNV  for  Pahute  events.  The  turning  ray 
model  for  propagation  is  assumed.  The  Mueller-Murphy  source  scaling  relations  and  a  t’  value  of 
0.1  s  are  also  assumed.  As  indicated  by  arrows,  the  observed  pP  arrival  is  late  but  approximately 
the  same  in  size  as  the  elastic  prediction. 
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Figure  30.  An  enlarged  map  showing  the  locations  of  stations.  Below  is  a  seismic  section  of  P 
waves  observed  from  the  Yucca  event  FLASK.  The  onsets  of  the  waveforms  are  relatively 
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event  FLASK.  The  onset  of  the  wave  form  appears  relatively  stable,  and  there 
is  clear  Pn  tnoveout  which  is  sufficient  reason  to  attempt  development  of  a  crustal 
mode 1 . 

We  began  our  numerical  experiments  by  attempting  to  use  a  layer  over  a  half 
space  model.  Synthetics  were  computed  with  generalized  ray  theory.  The  travel 
times  to  Las  Vegas  constrained  the  model  to  be  relatively  slow.  The  speed  of 
the  mantle  layer  appeared  to  be  only  7.9  km.  If  the  travel  time  was  matched  by 
slowing  the  crustal  layer,  the  contrast  across  the  Moho  became  too  large  and  the 
ratio  of  the  head  wave  to  the  reflected  energy  was  predicted  to  be  too  small. 
Introducing  a  gradient  into  the  crust-mantle  transition  appears  to  be  one  of  the 
best  among  several  possible  explanations.  A  transition  zone  of  about  7  km.  an 
average  crustal  speed  of  about  5.5  km/s  and  a  mantle  speed  of  7.9  km/s  appeared 
to  fit  the  observations  the  best.  Our  preferred  model  is  plotted  in  Figure  31. 

Unfortunately,  though  the  NTS -Las  Vegas  model  fits  the  relative  amplitude 
and  timing  of  the  data  in  Figure  30,  it  does  not  fit  the  later  complexity.  This 
was  disappointing  since  after  the  initial  3  or  4  strong  swings  the  later  coda 
is  relatively  weak  until  the  arrival  of  S.  At  this  juncture,  we  decided  to 
again  draw  on  our  experience  from  modeling  teleseismic  P  to  improve  our  physical 
understanding  of  Pn.  The  late  arrival  of  pPn  (Figure  28)  and  the  additional 
complexity  in  the  compressional  waveform  corresponded  exactly  to  features  we  had 
noted  previously  in  teleseismic  P.  We  had  developed  a  theory  and  performed 
modeling  experiments  to  show  that  these  features  could  be  explained  relatively 
simply  through  the  phenomenon  of  spall.  We  developed  this  model  in  working  with 
the  two  Amchitka  tests  MILROW  and  CANNIKIN  (Burdick  et  al.,  1984). 

Figure  32  shows  the  source  representation  we  developed  to  represent  the 
spalling  process.  We  employed  many  of  the  ideas  used  in  standard  dislocation 
representations  of  earthquakes.  As  indicated  in  the  Figure,  the  explosion  itself 
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Figure  32.  A  scale  drawing  of  the  models  for  Che  explosion  and  spalling 
process.  The  explosion  is  itself  an  isotropic  point  source  and  the  spall  is 
represented  by  a  shallow  disk  of  tension  dislocation. 
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is  represented  as  an  isotropic  point  force  with  an  appropriate  RDP  as  usual. 
We  then  represent  a  spall  source  as  a  shallow  disk  of  tension  dislocations.  The 
dislocations  are  effectively  replaced  with  a  distribution  of  force  systems  just 
as  an  earthquake  fault  is  replaced  with  a  distribution  of  double  couples.  Tension 
dislocations  can  be  expressed  as  a  weighted  sum  of  the  basic  double  couples  with 
the  weight  depending  on  the  elastic  constants  of  the  medium.  For  the  Amchitka 
tests,  strong  motion  data  were  available  from  the  spall  zone  from  which  we  were 
able  to  specify  the  kinematics  of  motion  on  the  tension  dislocations.  With  this 
simple  model  and  the  relatively  well  constrained  kinematics  of  motion  we  were 
able  to  explain  the  anomalies  in  the  MILROW  and  CANNIKIN  data. 

Figure  33  illustrates  in  a  simple  way  how  the  spall  model  can  account  for 
the  observed  delays  in  pP.  As  is  well  known,  the  spall  crack  opens  when  the 
tension  from  the  downgoing  surface  reflected  energy  exceeds  the  compressions 
generated  by  the  upgoing  wave  and  the  overburden  pressure  as  well  as  the  tensile 
strength  of  the  containment  medium.  As  this  disk  of  tension  dislocations  begins 
to  open,  it  radiates  downward  compressional  energy  coincident  with  the  downgoing 
tension  of  pP.  This  is  depicted  in  the  first  two  rows  of  Figure  33.  The  pP  and 
spall  opening  events  will  tend  to  cancel  each  other  to  some  extent  depending  on 
the  relative  strength  of  the  explosion  and  spall  sources.  The  spall  event  removes 
energy  from  pP  as  one  should  expect.  It  also  radiates  compressional  energy 
upward  which  reflects  back  into  the  earth  as  tension.  This  energy  arrives  after 
pP  from  the  explosion  with  the  same  sign.  It  was  our  suggestion  that  this  second 
spall  phase  was  the  phase  that  was  being  identified  as  a  late  arriving  pP.  As 
shown  in  the  figure,  there  will  also  be  a  set  of  arrivals  coming  later  associated 
with  the  spall  closing  event.  These  phases  cause  some  complexities  in  the 
waveform  that  are  in  fact  observed  in  the  Amchitka  data.  The  bottom  trace  in 
the  figure  shows  the  composite  of  the  radiated  pulses.  Of  course,  the  model  is 
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very  simplistic,  but  it  does  provide  a  framework  in  which  we  can  at  least  begin 
to  address  the  problem  of  the  effects  of  spall.  In  the  case  of  Amchitka,  the 
free  surface  records  indicated  that  the  sizes  of  the  explosion  and  spall  sources 
were  comparable,  so  this  is  an  issue  that  needs  to  be  addressed. 

The  results  of  our  waveform  modeling  studies  from  which  we  developed  the 
NTS  to  Las  Vegas  model  in  Figure  31  and  of  our  spall  modeling  study  are  shown 
in  Figure  34.  Four  representative  P  waveforms  from  Figure  30  are  shown  on  the 
right.  The  Pn  is  small  and  emergent.  The  phase  we  have  labeled  as  PmP  is 
composed  of  several  generalized  rays  just  past  critical  angle  in  the  lid  gradient. 
The  phase  labeled  pPmP  is  composed  of  the  associated  free  surface  reflected  rays 
and  they  are  of  course  opposite  in  sign.  The  first  row  of  synthetics  fits  the 
observed  waveforms  for  the  first  few  oscillations,  but  there  are  consistently 
at  least  two  continuing  strong  oscillations  in  the  data  which  are  not  predicted. 
The  center  row  of  synthetics  shows  the  effect  of  adding  in  the  spall  source. 
We  estimated  the  delay  of  pP  by  spall  at  0.3  seconds  based  on  measurements  from 
the  deconvolved  waveforms  like  those  in  Figure  28.  The  interference  of  the 
various  arrivals  was  assumed  to  be  as  illustrated  in  Figure  33.  The  amplitude 
of  the  spall  source  was  assumed  to  be  equal  to  the  explosion  based  on  our  Amchitka 
experience.  The  dwell  time  of  the  spall  was  adjusted  to  fit  the  data  with  a 
value  of  0.3  seconds  appearing  to  fit  best.  The  match  with  the  observations  is 
good  enough  to  suggest  that  spall  is  an  important  part  of  the  explosion  source 
which  should  be  taken  into  account.  The  fact  that  a  gradient  model  is  required 
suggests  that  the  P  wave  propagation  to  regional  distances  in  the  western  US. 
can  be  best  modeled  with  a  continuous  positive  gradient  (strong  near  the  free 
surface)  from  the  free  surface  down  to  the  lid  of  the  low  velocity  zone. 
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Figure  34.  A  comparison  between  the  data  and  the  synthetic  seismograms.  The  synthetics  on 
the  left  include  the  effects  of  the  structure  and  an  explosions  source.  Those  in  the  center 
include  the  effects  of  spall.  The  data  and  the  synthetics  are  in  good  agreement. 


CONCLDSIONS  AND  RECOMMENDATIONS 

In  overview,  the  most  surprising  and  important  result  o £  this  effort  is 
that  there  is  a  great  deal  more  organization  and  deterministic  character  in  high 
frequency  regional  data  than  one  might  have  expected  from  the  review  of  previous 
work.  There  is  no  doubt  that  there  is  much  scattering  resulting  in  unstable 
high  amplitude  coda,  but  if  one  focuses  on  specific  times  in  the  records  when 
wave  propagation  is  simple  to  understand  one  can  perform  meaningful  waveform 
modeling  in  the  time  domain.  In  this  investigation  we  have  worked  with  Pn  and 
the  onset  of  critical  reflection  from  the  Moho ,  but  it  seems  likely  that  related 
approaches  will  work  for  tangential  Sn  and  critical  SH  from  the  Moho. 

Other  more  specific  and  new  results  have  emerged  from  our  efforts  in 
developing  numerical  methods  and  performing  numerical  tests  with  them.  We  have 
developed  wavenumber  integration  as  a  practical  means  of  computing  Pnl  synthetics, 
though  it  requires  judicious  use  of  a  Q  model,  damping  through  integration  in 
the  complex  frequency  plane  and  wavenumber  filtering.  We  found  that  the  results 
of  the  Pnl  calculations  do  not  depend  strongly  on  the  Q  model,  but  it  is  numerically 
convenient  to  use  one  because  it  moves  poles  from  the  integration  contour  in  the 
k  plane.  Our  numerical  tests  showed  conclusively  that  an  appropriate  model  for 
the  velocity  gradients  near  the  free  surface  is  very  significant  to  the  accurate 
prediction  of  short  period  Pnl .  Conversely,  gradients  near  the  crust  model 
transition  are  not  important  to  waveshape  once  the  energy  trapped  there  is  well 
past  critical.  The  fact  that  one  simple  layered  model  (our  model  U)  did  a  more 
than  adequate  job  of  explaining  all  of  the  observations  from  Pahute  and  Yucca 
is  a  strong  indication  that  short  period  Pnl  is  a  very  stable  and  deterministic 
type  of  regional  energy.  This  is  promising  for  its  future  utility  in  solving 
the  discrimination  problem. 
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The  progress  in  understanding  the  short  period  Pn  waveform  over  the  course 
of  this  investigation  has  been  dramatic.  Not  only  does  it  appear  that  the  phase 
may  provide  the  core  of  a  physically  based,  transportable  discriminant,  but  it 
also  provides  a  new  seismic  window  through  which  to  observe  the  complex  source 
process  of  nuclear  event,  spallation  and  tectonic  release.  These  processes  have 
been  studied  intensively  over  the  last  decades  using  teleseismic  short  period 
P.  These  investigations  always  seemed  to  produce  inconclusive  results  because 
of  trade  offs  between  source  processes  and  Q.  It  appears  that  almost  nothing 
can  be  conclusively  resolved  if  the  possibility  of  laterally  varying  and  strongly 
frequency  dependent  Q  is  considered.  The  average  level  of  t*  appears  to  be  of 
Che  order  of  0.8  to  1.2  sec.  which  amounts  to  a  very  strong  effect  on  the  waveform. 
The  level  of  t*  for  Pn  appears  to  be  only  about  0.1s;  a  drop  of  an  order  of 
magnitude  in  the  attenuation  exponential.  This  affords  very  definitive  studies 
of  the  source  processes.  We  strongly  recommend  further  work  on  the  short  period 
Pn  waveform  discriminant  in  the  future. 


This  relatively  successful  two  year  study  could  not  have  been  carried  out  without 
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a  number  of  individuals  at  several  institutions.  In  particular,  we  would  like 
to  thank  Dr.  Steven  Taylor  at  LLNL,  Dr.  Hiroo  Kanamori  at  Caltech  and  Dr.  Duncan 
Agnew  at  Scripps  Institute. 
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